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ABSTRACT 


Landing aircraft on board carriers is a most delicate phase of flight operations at 
sea. The ability to predict the aircraft carrier’s motion over an interval of several sec- 
onds within reasonable error bounds may allow improvement in touchdown dispersion 
and a more certain value for a ramp clearance due to a smoother aircraft trajectory. 
Also, improved information to the Landing Signal Officer should decrease the number 
of waveoffs. 

This work indicates and shows graphically that, based on the data for pitch, heave 
and roll measured for various ships and sea conditons, the motion can be predicted well. 
The predictor was designed on the basis of Kalman’s optimum filtering theory for the 
discrete time case, adapted for real-time digital computer operation. 
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I. INTRODUCTION 


The landing phase of VIOL, STOL and HELO aircraft aboard a ship without a 
ship's Inertial Navigational System represents a complex operation and a demanding 
task. The last 10 to 15 seconds before the aircraft touchdown involves terminal guidance 
and control problems, where not only the aircraft is disturbed by several kinds of 
stochastic disturbances (wind), but also the touchdown point (on the ship) is being 
moved randomly. Despite the wind disturbances and the final point (target) random 
motion, the landing accuracy specified for ship operations is very high, i.e., a few tens 
of feet longtitudinal landing dispersion. Such a terminal point problem is made tractable 
in a most natural wav by assuming that the ship’s position can be predicted for several 
seconds ahead so that the aircraft is guided toward the future position of the touchdown 
point. The scope of this study was to establish to what extent a stochastic process, like 
the ship’s motion, is predictable over moderate periods of time. 

Graphical results obtained throughout this estimator’s feasibility study concerning 
the relationship between the ship’s states: heave, pitch and roll versus the ship’s esti- 
mated states; heave, pitch and roll; and the influence of noise measurement are pre- 
sented. Digital simulations show that the prediction accuracy was very accurate for 
different swell wave height and period. The feasibility of estimating the ship’s motion 
of acceptable bounds of error can also lead to improvement of the Landing Signal Of}i- 
cer’s (LSO) decision policy for waveoffs. 

This study was based on the use of the ship’s motion equations and the simulated 
measuring instrumentations existing on board ship in order to get the predicted motions. 
Using this information, a predictor based on Kalman’s theory of optimum estimation 
was designed. 

Several circumstances contributed to the success of this approach. The size and 
mass of the ship significantly filter the motion of the sea. A complete landing operation 
is short enough that the stochastic processes are reasonably taken to be stationary. 
Finally, the prediction interval is onlv a small fraction of the time it takes each aircraft 
for final approach. 

This work is divided into three parts: 1) the derivation of the mathematical model 
of the ship’s motion, 2) the rationale in the implementation of Kalman filter and pre- 


dictor equations, and 3) discussion of the simulated computer results obtained. Since 


we are interested only in the most critical aspects of the landing operation, namely the 
characteristics of the longtitudinal channel, we merely investigate in the sequel the pre- 
dictability of the heave, pitch and roll motion of the ship where yaw is omitted. 


Il. PROBLEM STATEMENT 


A. GENERAL 

The ship’s attitude estimation used in this thesis involves calculation of heave, pitch 
and roll of the ship moving in the general direction towards a swell while recovering an 
aircraft. The motion of the ship is given inz , @ and @ coordinates as shown in Figure 
e 

This problem will be developed using state space methods. Given the heave, pitch 
and roll motion (the measurements) received through the ship’s sensor system, we are 
interested in estimating the heave, pitch and roll of the ship. The state variables for this 
plant are z,,z,,0,, 6, , @, and Dy 


B. SHIP MOTION 

A ship moving on the surface of the sea is almost always in oscillatory motion. The 
different kinds of oscillatory motions that a ship experiences can be described with the 
help of Figure 1, which shows the six kinds of motion, three linear and three rotational 
about the principal axes. Accordingly, 


The six kinds of ship’s motions are: 


l. a 
vig 


2 
3. c = heaving -- motion vertically up and down. 
4 


surging -- motion backwards in the direction of the ship travel. 


swaying -- athwartship motion of the ship. 


. d = rolling -- angular motion about the longitudinal axis. When the ship rolls, it 
lists alternately from starboard to port and then back to starboard. 


Lo) 


e = pitching -- angular motion about the transverse axis. When a ship pitches, it 
trims alternately by the bow and by the stern. 


6. f = yawing -- angular motion about the vertical axis. 


x-axis is longitudinal axis 
y-axis is transverse axis 
z-axis is vertical axis 





Figure 1. The x, y, and z-axes of a ship 


Onlv three kinds of motion, namely heaving, pitching and rolling (z,@,@ ) are 
purely oscillatory motions, since these motions act under a restoring force or moment 
when the ship is disturbed from its equilibrium position. These ship motions are the only 
motions being considered in this thesis. 

Although in reality a ship experiences all six kinds of motion simultaneously, only 
one motion of the three angular motions will be treated at a time in the following 
sections. However, one must bear in mind that any one kind of motion is not inde- 
pendent of the others; in the first approximation, coupling between the motions is neg- 


lected in order to simplify the problem. 


C. SYSTEM MODEL 
The svstem to be modeled in this problem is the motion of a surface ship at sea 
during an aircraft recovery phase. This is a linear, time-angle system that can be de- 


scribed in a discrete equation of ship motion. The state space equation is 
Xpgy = Dydy + UV yuy (2.1) 


where 
xX,= parameter to be estimated (state vector) heave, pitch and roll 


,= state transition matrix which describes how the states of the dynamic system are 
related. 


I,= system random uncorrelated coefficient matrix. 


u,= sequence of random uncorrelated inputs. 


From Equation (2.1) and the above assumptions, the state vector for heave, pitch and 


roll is 


(AD) 


G- O68 D DNN 


and the system state equation can be approximated for T < < 1 


z i Ce © OO Mis 0 0 0 
Z =—¢, 1—b, 0 O =OMmNONIE te 070 
6 0 0 0 oO fie Fi (Org F3 
Ay |k+1 = y |k Uy We. 
6 0 0 -g¢ 1-b 0 O 6 0 Tg, 0 
ob > 0 °0 +O uel: 0 0 0 
db C0 0 OSes 0 O Ts, 

where 


u,=€, cos we + plant noise 

¢,= swell wave amplitude 

®,= swell angular frequency 

c, = g,= heaving restoring force coefficient 
b,= heaving damping force coefficient 

c, = g, = pitching restoring moment coefficient 
b, = pitching damping moment coefficient 

c, = g,=rolling restoring moment coefficient 


b,= heaving damping moment coefficient. 


The system noise process for the ship motion estimation problem is a function of the 


random uncorrelated coefficient matrix, [, , and the random forcing function, u,. 


D. MEASUREMENT MODEL 
For a linear measurement process, the measurements are linearly related to the state 


variables and can be modeled using the linear measurement equation 
Zp = HyXy + Vy (2.4) 


where 
z,= set of measurements. 


H,= observation matrix that gives the noiseless relationship between the measure- 
ments and the state vector. 


X,= state vector 


y,= sequence of random uncorrelated measurement noise. 


In this ship motion estimation problem, the measurements are the heave, pitch and 
roll angles taken from the ship sensors. 


Zz=Cl OJat+y (2.5) 


The linear equation with the measurement data available can now be processed with 
the Kalman filter. The type of noise that affect the ship motion data are assumed to be 
white gaussian noise. This noise model was used in the computer simulations and shown 
in Figure 2, with zero mean and variance of one. 

The complete development of this niodel can be found in Bhattacharyya [Ref. 1: 
pp.35-180], Wah [Ref. 2: pp.31-47], Korvin- Kroukovsky [Ref. 3] and Blagoveshchensky 
[Ref. 4]. 


The corresponding block diagram of the system is shown in Figure 3. 


Random Noise vs Time 
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Figure 2. White Noise Model 
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Figure 3. 9 Ship Motion Model Block Diagram 


Ii]. KALMAN FILTER THEORY 


A. GENERAL 

Filtering refers to the process of estimating the state vector at the current time based 
upon all past measurements. An optimal filter concentrates on optimizing a specific 
performance measure used to determine the quality of the estimate. The Kalman filter 
is an optimal filter in a class of linear filters that minimize the mean square estimation 
error between the actual and desired output. In other words, the Kalman filter attempts 
to minimize the elements along the main diagonal of the state error covariance matrix. 
The filter is a recursive algorithm for processing discrete measurements or observations 
in an optimal manner, [Ref. 5: p. 101] and [Ref. 6]. It requires a priori knowledge of the 
state estimate (x,_,) and its error covariance (P,_,), and the current observation (z,). The 
Kalman filter is the proper algorithm to be used when both the system model and the 
measurement model are linear functions of the state variables and these models can be 


described by the equations 
App = Deke + VU yUy (3.1) 
Zp = Ayky + Xe (3.2) 


In this equation, ®, and T, are constant matrices and H, is a linear constant func- 
tion of the state variable x,, while u, is the random forcing function and y, is the meas- 
urement noise. The plant random forcing function and measurement noise are assumed 


to be uncorrelated (white Gaussian) with zero mean. That is 


El(k).w'()] = O'(K)Sy (3.3) 
and 
ETv(k).v'()] = BK) (3.4) 
where 
Dy ell, BS) on 
=0, k#j 


B. NOISE PROCESSES 
The calculation of the error covariance matrix and the filter gain matrix requires the 
covariance matrices for the uncorrelated noise process uw, and vy, For the measurement 


noise process, y,, the covariance matrix is 
7 
ELy,v; ] = Ry (3.6) 


where R, is defined as the state measurement noise covariance matrix. It is based on the 
sensor accuracy and accounts for unknown disturbances such as steps, white noise, or 
imperfections in the plant model. The variance of the white noise model used in the 
computer simulations was one. 

The state excitation matrix, Q,, used in the Kalman filter represents the system noise 
process and is a function of the system noise coefficient matrix, [,, and the random 
forcing function, u,. This matrix is given by 


O,= (1 ,07.0 0 (3.7) 


where I, is the same as in Equation (2.3). The Q, matrix allows for any random heave, 
pitch and roll ship motion well as inaccuracies in the system model. The magnitude of 
QO, has a direct bearing on the magnitude of the state error covariance matrix and it 
prevents the covariance matrix from becoming singular by ensuring some uncertainty in 


the state estimates. 


C. INITIALIZATION AND OPERATION 

In the ship motion estimation, Kalman filters are used to minimize the ship’s atti- 
tude errors. Prior to processing the measurement data, the filter must be initialized with 
an initial state estimate and an initial error covariance matrix. This initialization process 
is a Very important step in the filter operation and gross inaccuracies in this step may 
cause the filter to diverge. Divergence occurs when the calculated covariance errors be- 
come much smaller than the actual covariance errors. This causes the estimate values 
of the states to pull away from the actual value and is also used as a startup value for 
the recursive scheme. 

The basic operation of the filter is a relatively straightforward recursive process. 
Based on Kirk [Ref. 6 : pp.4-1, 4-62] and Gelb [Ref. 5], the discrete time version of 


equations used in the Ka!:...n filter are: 


Seer) = P deen) (3.8) 


10 


Poe-iy = Px? “eye + Oy (3.9) 


Gy = Puen AePoe-tyHe + Ry (3.10) 
Xue ee Leet) + Gy(Zy — HyXcqx-t)) (3.11) 
Poguy = 1 — Gy) Panry (3.12) 


where 
Xue) = projected ahead state estimate 
@, = state transition matrix given by Eq. (2.3) 
Pus) = projected ahead state error covariance matrix 
Q, = state excitation covariance matrix given by Eq. (3.4) 
G, = Kalman gain matrix 


R, = state measurement noise covariance matrix given by Eq. (3.3) 


= 
| 


linearized measurement matrix given by Eq. (2.5). 


From the given equations, we can start the filter processing operation. The a prion 
state estimate are calculated using the © matrix shown in Equation (2.3) where 7 1s the 
sampling time in seconds between the observed ship motions (heave, pitch and roll). 
{It is assumed ship motions are received simultaneouslv through the ship sensors.) 

The Kalman gain matrix serves to minimize the mean square estimation error and 
is an indication of how much emphasis or weight will be placed on the current observa- 
tion. If P,,-;)18 small, the Kalman gain matrix will also be small due to the finite value 
of R,. If the P,,,-1 is relatively large, the gain is approximately one. By rewriting the 
equation for the calculation of the state estimate Equation (3.11), as 


Sug = (1 — GH Xetn—1y + Gre (3.13) 


we can see how the Kalman gain matrix directly affects the weight placed on the current 
observation z,. A large gain, indicating a large a priori error covariance, will place more 
weight on the current observation as the filter tries to correct the states. A small gain, 
indicating a small error covariance, places less emphasis on the new observation. 


If the Kalman gain is expressed as 


G, = hfe ule (3.14) 


11 


it can be seen that the gain matrix is “proportional” to the uncertainty in the estimate 
Py 4-1» and “inversely proportional” to the measurement noise R,. For a large R, and a 
small P,,-); the measurement in Equation (3.2) is due mainly to noise and only small 
corrections should be made in the state estimate. However, if R, is small and P,.,-1) is 
large, the measurement contains considerable information about the errors in the esti- 
mates and therefore, a strong correction should be made to the state estimates [Ref. 5: 
pp. 127-8}. 

The computation process for prediction is divided into the following steps: 1) cal- 
culate the estimate using Equation (3.11) and 2) use Equation (3.8) for the desired pre- 
diction. 

The block diagram of the system with Kalman Filter is shown in Figure 4 and the 
system predictor block diagram is shown in Figure 5. 
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Figure 4. Ship Motion Model Diagram With Estimator 
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Figure 5. 9 Ship’s Motion Predictor Block Diagram 


IV. COMPUTER SIMULATIONS 


A. GENERAL 

A measured simulated ship motion (heave, pitch and roll) and the forcing function 
data were chosen based from the ship speed of 25 knots, heading 000 degrees true and 
the swell wave heading direction at 120 degrees true. The sea state for a swell wave 
height of 11 feet is between 5 and 6 based on Bhattacharyya [Ref. 1: p. 104] which can 
be described that large waves begins to form while white crests are extensive everywhere 
(probably some spray). The wind force (beaufort) is classified as level 6 which is char- 
acterized by a strong breeze with the velocity of 25 knots. 

These measured data were used in the Matlab computer program ([Ref. 7] and [Ref. 
8]) algorithm (as shown in Appendix A) to generate the graphical results to see the 
phvsical representation of the ship motion oscillations. These graphical representations 
of swell, heave, pitch and roll, which are the results from the above computer runs, have 


the following parameters listed in Table 1 and are shown in Figures 6, 7,8 and 9. 


Table 1. WAVE AND SHIP MOTION PARAMETERS, TS =0.1 SEC 


[Coefficient | Swell | Heave | Pitch | Roll] 
ami) | - | - 
aa | me | on |e 
[Restoring [| 0 | osm] oss | 0391 _| 


ealmenmeGain OFl2355 0.1191 0.2197 
0.0841 0.0756 OAT PA 
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Figure 6. Swell Wave Motion with a Height of 11 feet 
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Figure 7. Ship’s Heave Motion caused by 11 feet high of Swell 
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Figure 8. Ship’s Pitch Motion caused by 11 feet high of Swell 
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Figure 9. Ship’s Roll Motion caused by 11 feet high of Swell 


With the simulated observed ship motion data for heave, pitch and roll generated 
from the above computer runs, the Kalman filter was used to calculate the estimates by 
using the Matlab computer program algorithm generated as shown in Appendix B. The 
resulting estimated data were represented in graphical form to show the difference be- 
tween the observed and the estimated ship motions. These graphs, with estimates, are 
shown in Figures 10, 11 and 12. Computer runs were also conducted to show the 
graphical representations of the predicted ship motions heave, pitch and roll. A Matlab 
computer program is shown in Appendix C for this calculation. These graphical pred- 
ictions of ship motions are shown in Figures 13, 14 and 15. As we can see, the estimates 
and predictions illustrate a well-defined representation of the ship motion that can be 
used in the prediction of the ship motion in the future time which is to be relayed to the 
Landing Signal Officer for his decision for waveoffs and to the aircraft's pilot through 
the aircraft's TACAN during the recovery phase of flight operations. 
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Figure 10. Observed and Estimated Heave Motion caused by 11 feet high of Swell 
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Figure 11. Observed and Estimated Pitch Motion caused by 11 feet high of Swell 
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Figure 12. Observed and Estimated Roll Motion caused by 11 feet high of Swell 
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Figure 13. Observed and Predicted Heave Motion caused by 11 fect high of Swell 


ay 


Ships Observed Pitch and Predicted Pitch 


0.08 


Figure 14. 


t 
! 
t 
t 
' 
' 
uo] 
@ 
> 
Ma 
el 
7) 
a 
° 





Predicted 
20 30 40 50 60 70 80 90 100 
Time(sec) 


10 


0.06 
0.04 
0.02 


(ped) ajduy Yoyd 


Observed and Predicted Pitch Motion caused by 11 feet high of Swell 
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Figure 15. Observed and Predicted Roll Motion caused by 11 fect high of Swell 
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V. CONCLUSIONS 


The feasibility of estimating an aircraft carrier motion at sea by measuring the ship’s 
actual position and motion was investigated. The ship’s motion representation was 
generated based on ship’s motion mathematical model taken from [Ref. 1]. Subse- 
quently, a Kalman filter-estimator adapted for real-time computation on a digital com- 
puter, was generated. The results obtained show that the desired prediction time can 
be reached with reasonably acceptable errors. 

Being able to predict the ship’s motion can lead to an improvement of aircraft 
landing accuracy and safety. This can be accomplished, for instance, by generating ter- 
minal guidance (landing) laws making use of the future ship’s motion and position. 
Moreover, the possibility of prediction can improve the LSO information and policy for 
landing acceptance or wave-offs. The possibility of processing the measured motion by 
Fast Fourier Transform Algorithms (FFT), in order to obtain the estimated ship motion 
parameters such as the heave ({ z ), pitch (8) and roll (@), in real time may lead toward 


an adaptive predictor scheme. 
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APPENDIX A. MATLAB PROGRAM | 


This Matlab computer program calculates the Observed Ship’s Motion such as the 
heave, pitch and roll of the ship. Required coefficients must be available to run the 
program. 


% A Matlab program to generate the ship's heave, pitch and roll 
% roll motion with a given Swell wave force of N sin ws*t. 


Rrpbbbierrbioiioerwobiereooiee SWELL parameter hbo 010/110 1/0/10 1010010 
% 


% 
% 
diary p.m 
% Swell period (Ts) seconds 
Ts=7; 
% Swell angular frequency 


ws= 2*pi*1/Ts ; 


% Wave Height (N) feet 
N=11; 


> 


bi=N*ws; 


% Restoring coefficient 
cs=ws**2 : 

% 

ibepebipbepppoppomeoereeee Heave parameter %hhhpbrpwpmmm iors /0/e 0/07 

% 


% Heave period (th) seconds typical for carriers 
tn— 7 ls 

% 

% Heave angular frequency (wh) rad/sec 


wh= 2*pi*1/th : 


% 


% Heave restoring cefficient (cz) 
cz=wh 
% Heave damping factor (zetah) 
zetah= 0. 6; 


% Damping coefficient (bh) 
bz=2*zetah*wh § ; 
% 
% 
Lhhppepoprorhoppphporeeooe Pitch parameter %L%hwow oooh o/010 7010 1010101010 100 
id 


a 


2 


% 
% Pitch period (tp) seconds typical for carriers 


tp=10; 

% Pitch angular frequency (wp) rad/sec 
wp= 2*pi*l/tp g 

% 

% 

% Pitch Restoring coefficient 
cp=wp 

% 

% Heave damping factor (zetap) 

zetap= 0.6 : 
% 


% Damping coefficient (bp) 
bp=2*zetap*wp ; 


% 
% 
hhhhhhphhphhhhphhbhppphee ROlL parameter %hhbbphbbniniewoioie nies leteseh ble hele 
% 
to Roll period (tr) seconds typical for carriers 

tr= 18; 
% Roll angular frequency (wr) rad/sec 

wr= 2*pi*1l/tr 2 
% Roll Restoring coefficient 

cr=wr 

% 
% Roll damping factor (zetar) 


zetar= 0.6 ; 


% Damping coefficient (br) 
br=2*zetar*wr : 


% 
% 
Powbbwbbpeoewowovopoeseo Parameters of the system %%6%10 000% 0/00 10/0404 
% 


% where 

% =[];-- the A martix coefficient in continous time 
% b=[];-- the B matrix coefficient in continuos time 
% 

% Swell Parameters 


oe TSS SSS SSS SS SS SS 


as=-0 1;-cs 0-; 
bs=[ 0; 1]; 


% Heave Parameters 


gy =———. 
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ah=[(0 1; -ez**2 -bz] ; 


bh=[ 0; . 3]; 
% 
% Pitch Parameters 
4 ee 
% 
ap=(0) 1scp**2"=bpy 
bp=[0;. 3] ; 
% 
% 
qo Roll Parameters 
ve ———— 
% 


ar—(O 1: -cr**2Z =br] ; 
br=[0;. 6]; 
diary off 








coefficient of the system %%%%w»o 


ic 
% 
% Measurement Parameters 
Os 2] SSS SS SSS SS 
% 
c=[ 1 0); 
d=[ 0] ; 
dt=0. 1; 
q=eye(2); 
g=eye(2); 
r=1/10000; 
stat ab (2/5 
cs=[N 0]; 
% 
hhbbbbobppooos Determine the discrete 
% The Phi's and Gamma's of the system 
o 
% 


{ phih, gamh]=c2d(ah,bh,dt) ; 
[phip,gamp]=c2d(ap,bp,dt) ; 
{ phir, gamr] =c2d(ar,br,dt) 

[phis,gams]=c2d(as,bs,dt) ; 


* 
3 


wmpwpeeeeee Generate Input for the 
t=(0:.1:100)'; 
rand('normal'); 
v=sqrt(r)*rand(1000,1); 

% 

% 
us=dlsim(phis,gams,cs,d,v); 

% 

% 

hbebobeeeeeee Generate output values 

% 

% 
yh=dlsim(phih, gamh,c,d,us); 


system %%%io/0/0/0/0/010 10/010 16 0/010 f0 10 1010 10 0 1010 10 


Time axis 

Random are uniformly distributed 
in the interval (0.0,1.0). Random 
('normal') switches to normal with 
mean of 0 and variance of 1 


for the system %%%%%%0/070 1010 1010 10 1010 10.0 Jo "0 


co 


% 
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Lhhhhhhheeeore Plot the generated output vs time of the system 


te 


yp=d1sim(phip, gamp,c,d,us); 
yr=dlsim(phir, gamr,c,d,us); 


plot(yh), title(' Ships Heave Motion' jeg 
xlabel('Time (sec)'),ylabel('Heave angle (rad)'); 
delete pitch. met 

meta pitch 


plot(yp), title(' Ships Piteh Motion’ ae 
Xlabel('Time (sec)'),ylabel('Pitch angle (rad)'); 
meta 
plot(yr), eitiee Ships Roll Motion’ ); 
xlabel('Time (sec)'),ylabel('Roll angle (rad)'); 
meta 
plot(v) ,title(' Random Noise vs Time'); 
xlabel('Time (sec)'),ylabel('Noise Magnitude’ ); 
meta 
plot(us), Eicle( Swell with White Noise vs Time'); 
xlabel('Time (sec)'),ylabel('Swell angle (rad)'); 
meta 
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APPENDIX B. MATLAB PROGRAM 2 


This Matlab computer program calculates the Observed Ship’s Motion and esti- 
mates in real time such as the heave, pitch and roll of the ship. Required coefficients 
must be available to run the program. 


% This Matlab computer program calculates the Observed and Estimated 
% state of the Ship s Heave, Pitch and Roll based on the given data. 

% 
% 

Wiiieieriepienieprpppiwwemeeeee SWELL parameter %LLbioo0%0010/0070/070/0/0/010/0 10/070 
% 


% 


diary p.m 
% Swell period (Ts) seconds 
Ts=7 ; 
% Swell angular frequency 


ws= 2*pi*1/Ts  ; 


% Wave Height (N) feet 
N=11; 


% Restoring coefficient 
cS=ws**2 ; 
% 
% 
Db bpbhpprpppppbpppbpbeeee Heave parameter %%ALL%o% 1701 1970 70.10% 707010 %0.%0 10.10.10 
% 


% Heave period (th) seconds typical for carriers 
th= 11 2 

% 

% Heave angular frequency (wh) rad/sec 


wh= 2*pi*1/th : 
% 


% Heave restoring cefficient (cz) 
cz=wh : 
% Heave damping factor (zetah) 


zetah= 0.6; 


% Damping coefficient (bh) 
bz=2*zetah*wh ; 
% 
% 
LLM ppowpepmeooopeooes Pitch parameter  %%~K% 770101770. 707070. 10.70 10 10 Yo 10 10 10 


% Pitch period (tp) seconds typical for carriers 


at 


tp=10; 


%o Pitch angular frequency (wp) rad/sec 
wp= 2*pi*l/tp ; 

% 

% 


% Pitch Restoring coefficient 
cp=wp ; 
% 
% Heave damping factor (zetap) 
zetap= 0.6 : 


% 
% Damping coefficient (bp) 
bp=2*zet ap*wp ; 
we 
Wh 
Lhbbbhhabohbbbhhbbhhoohhn ROL parameter ME Gane aLue meus on canine ome 
Uh 


% 


% Roll period (tr) seconds typical for carriers 
tr= 18; 
% Roll angular frequency (wr) rad/sec 


wr= 2*pi*l/tr ; 


% Roll Restoring coefficient 
cr=wr ; 
% 
% Roll damping factor (zetar) 
zetar= 0.6 ; 


% Damping coefficient (br) 
br=2*zetar*wr ; 
% 
Lhbbbbhbhpphhbbbaphpoohoee Parameters of the system %ALLbbohobboobbhoeph 
% 


% where: 

% a=[];-- the A martix coefficient in continous time 
% b=[];-- the B matrix coefficient in continuos time 
ay 

/o 

% Swell Parameters 


of ——=_ = 











as=-0 1;-cs 0-; 
fs=[0; 1]; 


% Heave Parameters 


va SS oS SS SS 


ah=[0 1; -cz**2 -bz] ; 
£h=[/0;= 3]: 
% 


% Pitch Parameters 





% 

ap=[0 1; -cp**2 =bp] ; 

fp=[0;. 3]; 
% 
% 
% Roll Parameters 
G4 ee ee 
% 

ar—(O Ws=cr**2 =br] ; 

fr=[0;. 6]; 
diary off 
% 
% 
% Measurement Parameters 
v5 ee 
% 
% 

c=[1 0]; 

dt=0. 1; 
% 
%o Swell Noise 
ese ee ae ee 

q=1/100; 
% 
% Measurement Noise 
gS Sa ee 

r=1/8500; 

d=[ 0] ; 
%o 
% 
Dhwwmweoeeeooeeee Determine the Discrete coeficient of the system %%%% 
Yo 
% Swell 
y% 
% 
[ phis, gams] =c2d(as,fs,dt) p 
% 
% HEAVE 
 - 
% 
[ phih, gamh] =c2d( ah, fh,dt) : 

% 
% PITCH 
EES 
%o 
i [ phip, gamp] =c2d( ap, fp,dt) ; 
% ROLL 
% oe 


[ phir, gamr] =c2d(ar,fr,dt) : 


5) 


% 

Yo Fo 70:10 fo fe YoYo Yo Wo 

% 
kh=dlqe(phih, gamh,c,q,r) ; 
kp=dlqe(phip,gamp,c,q,r) ; 
kr=dlqe(phir,gamr,c,q,r) ; 


Calculate the Kalman Gains for Heave, Pitch and Roll %%%%%%% 


Lhpbophoppowrhpowhpeeoe Initial Condition ZAK hwmmm bohm —o% 


% 
xh=[ 0. 05; 0. 05] 
xhe=[ 0. 
xp=[ 0. 
xpe=[ 0.05; 0. 05] 
xr=[0. 
xre=[ 0.05; 0. 05] 
xs=[0.05;0.05] ; 
% 


hhrviohwopporwoeeoheee Random Noise Generators 2% brbmbbbmbh bio hle 


% 
rand('normal' ) 
w=sqrt(q)*rand(1,1500); 
v=sqrt(r)*rand(1,1500); 
t=0:.1: 100; 

% 

Yo 707010 70 10 1010 10.10.1010 10 70 10 10 70 10 10 fo 
for k=1: 1000; 

% 


Simulation Loop 2%: //o/0 /e 1001010 /0 1010/00 1010/0 


fowpperooioiopeooooooe iMmput for Heave, Pitch nad Roll %A%K%EE%% 


% 
xs(: ,k+1)=phis*xs(: ,k) + gams*w(k); 


Dbhphoepoe State calculation for Heave, Pitch and Roll %%%%%%%%% 


% 
xh(: ,k+1)=phih*xh(: ,k) + 
xp(: ,k+1)=phip*xp(: ,k) + 
xr(: ,k+1)=phir*xr(: ,k) + 


gamh*xs(k) ; 
gamp*xs(k) ; 
gamr*xs(k) ; 
% 


oO 


/o 
bhbhbehooo Output calculation 


% 
yh(: 
yp: 
yr: 


»kt+1)=c*xh(: ,k+1) + 
enlace p(s Kt) + 
ykt+t1)=c*¥xr(: ,k+1) + 


v(k+1); 
v(k+1); 
v(k+1); 
eo 
a 
Joie fo Jo10 10 10701 


° 


State estimate calculation for 


xhe(: ,k+1)= phih*xhe(: ,k) + kh*(yh(: 
xhe(: ,k) )); 
xpe(: ,k+1)= phip*xpe(: ,k) + kp*(Cyp(: 
xpe(: ,k) )); 
xre(: ,k+1)= phir*xre(: ,k) + kr*(yr(: 
xre(; ,k) )); 
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for Heave, Pitch and Roll in Discrete form 


Heave, Pitch and Roll %%%%%%%%%% 
p kt 1) =e*( pha... 
k+1)=c*(phip™, . . 


»Krl)-e%( phir... 


end 


% 
% 


~o 
eo 


plotGeexh( lt. xhec).: )).,titlet Ships 

Observed Heave and Heave Estimate’ ) 
xlabel('Time(sec)'),ylabel('Heave Angle (rad)') 

delete esth. met 

meta est 


Plleect-sp(i,: ietsxpell,: )),title('Ships 
Observed Pitch and Pitch Estimate’ ) 


xlabel('Time(sec)'),ylabel('Pitch Angle (rad)') 
meta 


DlOtteeertl,: ),t;,xre(1,:)), title(' Ships 
Observed Roll and Roll Estimate’ ) 


xlabel('Time(sec)'),ylabel('Roll Angle (rad)") 
meta 
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APPENDIX C. MATLAB PROGRAM 3 


This Matlab computer program calculates the Observed Ship’s Motion and predicts 
in real time such as the heave, pitch and roll of the ship. Required coefficients must be 
available to run the program. 


% This Matlab computer program calculates the Observed and Predicted 
% state of the Ship s Heave, Pitch and Roll based from the given data. 
% 

% 

to 
dbp bb phpprrbobpbiphhhbeee SWELL parameter %RLLLAMbwmwmwmmowiwmweerwrerorere rer 
% 

% 


% Swell period (Ts) seconds 
Ts=7 ; 
% Swell angular frequency 


ws= 2*pi*1/Ts  ; 


% Wave Height (N) feet 
N=11; 


% Restoring coefficient 
cS=ws**2 ; 

% 

70 

iebbeeppprrwebpbbppphwheaeees Heave parameter %LLLLbwpmmbimmibprppeip rer ler 

% 


% Heave period (th) seconds typical for carriers 
th= 11 : 

% 

% Heave angular frequency (wh) rad/sec 


wh= 2*pi*1/th : 


% Heave restoring cefficient (cz) 
cz=wh f 
% Heave damping factor (zetah) 
zetah= 0. 6; 


% Damping coefficient (bh) 
bz=2*zetah*wh ; 
% 
% 
Lhpbppiewpweppweppeebpeees Pitch parameter  %~%W% Ibo /o1o 10/070 10.10 /0 10.10 701010 10 10.16 fo 


% Pitch period (tp) seconds typical for carriers 
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tp=10; 


% Pitch angular frequency (wp) rad/sec 
Wp= 2*pi*l/tp Q 


% Pitch Restoring coefficient 

» cp=wp ; 

% Heave damping factor (zetap) 
zetap= 0.6 g 

% 


te Damping coefficient (bp) 
bp=2*zetap*wp ; 


LEUUUALULUULALLALUGS RO1L parameter  GUU%LURGLLUUIUIALUAULA ALLS 
% 
% 


% Roll period (tr) seconds typical for carriers 
tr= 18; 
% Roll angular frequency (wr) rad/sec 


wr= 2*pi*l/tr ; 


% Roll Restoring coefficient 
cr=wr ‘ 
% 
vs Roll damping factor (zetar) 
zetar= 0.6 ; 


% Damping coefficient (br) 
br=2*zetar*wr ; 
% 
% 
Dhiohpeeepppppepphoeeorooee Parameters of the system %%LL%bwbpmwioieioilalaio 
% 


% where: 

‘% a=[];-- the A martix coefficient in continous time 
% b=[];-- the B matrix coefficient in continuos time 
9, 

% 

% Swell Parameters 





as=[0 1;-cs 0]; 
fs=[0; 1]; 


% Heave Parameters 


oy =p ee ee eee ee 


ah=[0 1; -cz**2 -bz] ; 
fh=[ 0; . 3]; 
% 
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% Pitch Parameters 


ap=[0 1;-cp**2 -bp] ; 
fp=(0;. 3]; 


% Roll Parameters 


WA — 





ar | ete cr 2 = ba |e 
fr=[0;. 6]; 





% Measurement Parameters 
Vb, SS SSSs 2 SS SSS SS SS SSS See 
% 
% 
el OD 
dt=0. 1; 
% 
% Swell Noise 
y «- Sencecooeess 
q=1/100; 
% 
% Measurement Noise 
yY ~) Scene eceee cic ooee 
r=1/3500; 
d=[ 0]; 
% 
% 
Lhbbbhppphhhabopoe Determine the Discrete coeficient of the system %2%% 
% 
% Swell 
ee were 
% 
[ phis, gams] =c2d(as,fs,dt) : 
% 
% HEAVE 
5 |) eee 
% 
[ phih, gamh] =c2d(ah, fh, dt) ; 
% 
% PITCH 
ye ee 
% 
[phip, gamp] =c2d( ap, fp, dt) 
% 
% ROLL 
Sh, enocs 
% 


[ phir, gamr] =c2d(ar, fr,dt) : 
% 
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% 
hbbpeeeeee Calculate the Kalman Gains for Heave, Pitch and Roll %%%%%%% 
% 
kh=dlqe(phih, gamh,c,q,r) ; 
kp=dlqe(phip,gamp,c,q,r) ; 
kr=dlqe(phir,gamr,c,q,r) ; 
% 
hhbhhhhhhoopoppphphhhr Initial Condition %%%%%%%%%%%%%%% moo % 
% 
xh=[-0.05;-0.05] ; 
xhe=[ -0.05; -0. 05] ; 
xp=[ 0.05; -0.05] ; 
xpe=[-0.05;-0.05] ; 
xr=[ 0.050; -0.050] ; 
xre=[-0.05;-0.05]_ ; 
<salOs05;)-0.05]: ; 


hhprbbwpeovivrieeoovoooee Random Noise Generators %rrboppmmiowb poole 
% 


rand('normal' ) 

w=sqrt(q)*rand(1,1200); 

v=sqrt(r)*rand(1,1200); 

(e=0s . ile Se 
% 
hhh phhhhhppowpebprhoe Simulation Loop ~~ rororo oo lore rorororo we lorore reo 
% 

for k=1: 950; 
% 
Lhbphhwhepopopeeeeoe Input for Heave, Pitch nad Roll %2ARK%%%% 
% 

xs(: ,k+1)=phis*xs(: ,k) + gams*w(k); 


or 


Mibeheene state calculation for Heave, Pitch and Koll %%X%%%%%%% 


xh(: ,k+1)=phih*xh(: ,k) + gamh*xs(k) ; 
xp(: ,k+1)=phip*xp(: ,k) + gamp*xs(k) ; 
xr(: ,k+1)=phir*xr(: ,k) + gamr*xs(k) ; 


% 
Peeoee Output calculation for Heave, Pitch and Roll in Discrete form 
% 
yh(: ,k+1)=c*xh(: ,k+1) + vCk+1); 
Veo: kl J—c*xp(: ,k+1) + v(k+1); 
vec. kel j=e*xr(: ,£kt1l) + v(k+1); 
% 
LbEhhbom State estimate calculation for Heave, Pitch and Roll %nenKnrnn% 
% 
xhe(: ,k+1)= phih*xhe(: ,k) + kh*(yh(: ,k+1)-c*(phih*... 
xhe(: ,k) )); 
RpeC et) Serehiptxpe(: 5k) + kp*(ypC: ,kt1)-c*(phip*... 
xpe(: ,k) )); 
xre(: ,k+1)= phir*xre(: ,k) + kr*¥(yr(: ,k+1)-c*(phir*®... 
xre(:,k) J); 
% 
% 
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% 

LE~%%% State Prediction Calculation for Heave, Pitch and Roll %%%Q%%%% 

% 
xhp(: ,k+1)=phih*xhe(: ,k) + gamh*xs(k) ; 
xpp(: ,k+1)=phip*xpe(: ,k) + gamp*xs(k) ; 
xrp(: ,k+1)=phir*xre(: ,k) + gamr*xs(k) ; 

end 
% 
% 


plot(t,xh(1,:), - ,t.xhp(l,: )), title Shins 
Observed Heave and Predicted Heave') 
xlabel('Time(sec)'),ylabel('Heave Angle (rad)') 
delete predicth. met 
meta predicth 


plot(t,xpG1;:),'-',t.xpp( le: )dotitle( Ships 

Observed Pitch and Predicted Pitch') 
xlabel('Time(sec)'),ylabel('Pitch Angle (rad)') 
meta 


plot(t,xr(l,:7),. -',t,xrpll.: Je tielec Ships 

Observed Roll and Predicted Roll’) 
xlabel('Time(sec)'),ylabel('Roll Angle (rad)') 
meta 
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